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ABSTRACT 

The  diffusion  and  recombination  process  in  an  ensemble  of  isolated  single  pairs  of  opposite  charges  is  usually 
described  by  the  Debye  - Smoluchowski  equation.  The  present  work  is  an  overview  of  a series  of  computer  simulations  of 
diffusion  and  recombination  of  ions  in  solution  performed  with  the  aim  to  determine  the  limits  of  validity  of  the  Debye  - 
Smoluchowski  equation. 

In  the  first  part  of  the  project,  the  calculations  were  performed  for  the  media  with  the  short  mean  free  path  (MFP)  of  the 
free  movement  of  ions  between  scattering  events,  i.e.,  for  the  conditions  of  the  diffusion  model  of  the  ion  transport.  Results 
were  obtained  on  the  probability  of  ion  survival  as  a function  of  time  and  the  probability  of  ion  escape  from  recombination 
at  infinite  time.  The  recombination  processes  in  the  clusters  of  non-separable  ion  pairs  and  the  bulk  recombination  of  ions  in 
solution  were  simulated.  The  deviations  of  the  multi-pair  kinetics  and  escape  probability  from  the  corresponding  results  of 
the  calculations  performed  on  a basis  of  the  Debye-Smoluchowski  theory  are  significant  but  we  found  that  the  Debye- 
Smoluchowski  recombination  rate  constant  can  be  applied  for  all  concentrations  of  ions. 

We  also  consider  the  effects  of  restricted  geometry  and  anisotropy  of  the  medium  on  the  kinetics  of  the  recombination 
of  oppositely  charged  species  in  each  other's  field  and  the  escape  probability.  This  model  roughly  corresponds  to  the 
electron  - cation  recombination  in  organic  crystals. 

The  results  for  short  MFP’s  were  then  compared  with  the  results  of  the  calculations  of  the  electron  - cation 
recombination  in  the  systems  characterized  by  long  MFP  of  electron  between  the  scattering  events.  It  was  found  that  for 
long  mean  free  paths,  above  10-20  % of  Onsager  distance  in  the  medium  {Rc=e2/4n££okT,  eeo  is  the  dielectric  permittivity  of 
the  medium,  e - electron  charge,  kBT  - the  thermal  energy),  the  recombination  theory  based  on  the  Debye  - Smoluchowski 
equation  is  no  longer  valid,  even  for  isolated  ion  pairs. 

Keywords:  diffusion,  ion  recombination,  rate  constant,  ionic  solutions,  computer  simulation. 


1.  INTRODUCTION 

The  present  theoretical  basis  for  diffusion-controlled  recombination  of  ions  is  provided  by  the  Debye  - Smoluchowski 
(DS)  equation  x 3 (see  1 for  a source  of  references) 

dC/dt  = V{D[  VC  + (C/kBT)  VUJ}  (1) 

where  C(r,t)  denotes  the  density  of  ions  as  a function  of  their  spatial  coordinates  and  time,  V is  the  interaction  potential 
between  ions,  D is  the  sum  of  the  diffusion  coefficients  of  cations  and  anions,  kB  and  T - the  Boltzmann  constant  and 
temperature,  respectively. 
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Unfortunately,  this  equation  describes  in  feet  only  a single  pair  of  interacting  particles  in  space.  If  the  interaction 
potential  V corresponds  to  the  Coulomb  interaction  -e7/4nrs£0  in  a dielectric  continuum  with  the  permittivity  ££m  then  we 
can  write  the  specific  version  of  the  equation  for  ions.  In  radial  coordinates  we  obtain 

dC/dt  -D0C  /dr*  + (2/r  + RJr2)  dC/dr } (2) 

with  the  boundary  conditions  at  the  encounter  radius  if  in  the  form  of  either 

C(Rt)=  0 or  (3) 

kRC(R,t)  = 4nR2D  (dC/dr  ) | ,=R  (4) 

where  the  former  is  known  as  the  Smoluchowski  condition  and  the  latter  as  the  radiation  boundary  condition,  with  kR  the 
reaction  rate  at  encounter.  The  parameter  Rc  is  called  the  Onsager  radius  and  is  defined  as  the  distance  at  which  the  energy 
of  the  Coulomb  interaction  in  dielectric  continuum  becomes  equal  to  the  thermal  energy  kBT : 

Rc=e?/4n££jcBT  (5) 

The  DS  equation  looks  deceptively  simple.  Its  approximate  solutions  are  known  already  quite  long  time  (see  1 for  a 
review)  but  the  fell  state-of-art  solution  for  C(r,t)  has  been  published  by  Hong  and  Noolandi  4 . When  the  spatial  variable 
from  C(r,t)  is  integrated  out,  we  obtain  P(t),  the  survival  probability  of  an  ion  pair  as  a function  of  time.  Proper  averaging  of 
P(t)  over  the  initial  spatial  arrangements  of  charged  species  leads  to  a time  dependent  quantity  that  can  be  compared  with 
the  macroscopic  recombination  kinetics.  Hie  limiting  value  of  the  spatial  integral  of  C(r,t),  at  infinite  time,  represents  the 
probability  that  an  ion  escapes  recombination  and  is  known  as  the  ion  escape  probability.  The  (simplified)  formula  given  by 
Onsager  1,4  is  as  follows 


Pesc  = {exp(-RJRJ  - exp(-R/R)}  /{l-  exp(-R/R)}  « exp(-R/RJ  (6) 

where  R0  is  the  initial  distance  between  the  counterions. 

The  DS  equation  has  been  derived  for  charged  species  with  a short  mean  free  path  (MFP)  between  consecutive 
collisions  with  the  molecules  of  an  isotropic  medium,  i.e.  for  ionic  species  or  for  excess  electrons  with  a relatively  low 
mobility.  Therefore,  incorrect  results  may  be  obtained  when  the  DS  equation  is  applied  to  a physical  system  for  which  one 
or  more  of  these  assumptions  is  violated.  In  this  papa"  we  consider  the  following  major  causes  of  breakdown  of  the  DS 
equation: 

(1)  Presence  of  multiple  ion  pairs  (which  leads  to  breakdown  of  the  scaling  properties  of  the  DS  equation); 

(2)  Anisotropy  of  diffusion  tensor 

(3)  Long  mean  free  paths  of  mobile  charges  and  velocity-dependent  scattering  mechanism  of  mobile  charges 

The  present  work  is  an  overview  of  a series  of  computer  simulations  of  diffusion  and  recombination  of  ions  in  the 
systems  containing  isolated  ion  pairs  or  a number  of  nonseparable  ion  pairs,  the  systems  with  short  or  long  MFP,  and  the 
systems  with  significant  anisotropy  of  diffusion  coefficients. 


2.  COMPUTER  SIMULATION  METHOD 

The  method  which  has  been  selected  for  calculations  of  the  diffusion  processes  belongs  to  the  class  of  the  Stochastic 
Dynamics  methods  (Allen  and  Tildesley  12).  The  method  ignores  the  neutral  particles  of  the  medium  and  treats  them  as  a 
dielectric  continuum  - only  the  ions  in  the  sample  are  treated  on  the  molecular  level.  The  influence  of  the  neutral  molecules 
is  simulated  by  the  random  walk  of  the  ions.  This  method  has  been  successfully  used  to  simulate  the  ion  recombination  in 
the  ionization  tracks  of  high-energy  particles  in  nonpolar  media  6 ■ 7’  13  . This  section  will  be  devoted  to  a more  detailed 
description  of  the  computer  simulation  programs. 

2.1.  The  diffusion  type  of  transport  of  charged  species  (short  MFP) 

The  initial  positions  of  the  ions  are  represented  by  a set  of  numbers  in  the  computer  memory.  Different  rules  of  the  ion 
distribution  in  the  systems  can  be  applied  - from  felly  random  distribution  of  ions,  through  simulation  of  the  Boltzmann 
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distribution,  to  regular  geometrical  arrangements  of  the  ions.  The  simulation  method  does  not  account  for  the  actual  velocity 
of  the  ions  nor  the  movement  of  the  neutral  molecules  of  the  medium.  We  assume  that  the  frequency  of  collisions  of  an  ion 
with  the  medium  molecules  is  so  high  and  its  mean  free  path  between  the  collisions  is  so  short  that  the  ion  velocity  is  a 
random  variable.  The  influence  of  the  medium  molecules  on  a given  ion  is  modelled  by  the  random  walk  (Brownian 
movement)  of  an  ion  in  continuous  medium.  The  random  vectors  are  obtained  from  the  random  number  generator  and  they 
are  properly  scaled  using  the  Einstein  relation  - the  mean  value  of  the  squared  distance  in  the  Brownian  motion  is 
proportional  to  the  time  of  motion  and  the  proportionality  coefficient  is  calculated  from  the  experimental  diffusion 
coefficient  <(Arf>  = 6DAt . The  ion  drift  vectors  in  electric  field  are  obtained  by  solving  the  electrodynamic  equations  of 
motion  using  experimental  mobilities  of  the  ions.  It  is  of  course  possible  to  include  any  type  of  external  field  affecting  the 
movement  of  ions. 

2.2.  The  transport  of  charged  species  with  long  MFP 

The  initial  conditions  in  the  long  MFP  case  include  also  the  velocities  of  the  electrons  which  are  sampled  from  the 
Maxwell  - Boltzmann  distribution  for  an  assumed  temperature.  The  random  directions  of  the  electron  velocities  are 
assumed.  The  simulation  method  does  not  account  for  the  actual  velocity  of  the  positive  ions  nor  the  movement  of  the 
neutral  molecules  of  the  medium. 

New  positions  of  the  cations  and  electrons  are  calculated  in  discrete  time  steps  using  the  two-timescale  technique.  The 
long  time  step  Atc  is  selected  for  cations.  Due  to  high  electron  mobility,  the  position  of  an  electron  can  change  quite 
significantly  during  the  cation  time  step.  The  short  time  step  Ate  is  selected  for  a given  electron  in  the  following  procedure. 
First,  an  actual  value  of  the  electron  path  X is  sampled  from  the  exponential  distribution  of  X.  Having  selected  the  time  step 
Ate,  the  trajectory  of  a given  electron  is  calculated  by  the  numerical  solution  of  the  electrodynamic  equations  of  motion  (for 
a classical  particle).  At  the  end  of  the  trajectory  the  mobile  charge  is  scattered.  The  scattering  mechanism  was  based  on  the 
assumption  of  complete  randomization  of  velocity  of  the  electron  on  collision.  The  model  seems  to  be  oversimplified, 
particularly  in  the  case  of  elastic  scattering  and  the  correlations.  On  the  other  hand,  the  assumption  of  the  full  randomization 
of  the  electron  momentum  makes  possible  the  direct  comparisons  between  the  diffusion  transport  and  the  long  MFP  cases 
and  systematic  investigation  of  the  recombination  kinetics  and  the  escape  probability  as  functions  of  the  MFP  length. 

For  the  both  cases,  short  and  long  MFP,  new  positions  of  the  ions  are  calculated  in  discrete  time  steps.  If  any  ion  pair 
fullfills  the  assumed  recombination  criteria  the  algorithm  simulates  the  ion  recombination  - the  ions  disappear  from  the  file 
of  the  "living"  ions  and  their  positions  are  no  longer  modified.  The  program  records  the  time  from  the  begining  of  the 
simulation  run  which  will  serve  to  construct  the  kinetics  of  the  recombination  process. 

The  simulation  runs  until  a prescribed  number  of  time  steps  (usually  104  or  105  time  steps)  or  is  terminated  in  the  case 
of  disappearance  of  all  the  ions  in  the  simulated  sample.  The  program  starts  again  for  a new  geometrical  configuration  of 
ions  and  is  repeated  for  a large  number  of  the  configurations.  To  obtain  a good  statistics  of  the  recorded  data  and  smooth 
kinetic  curves,  as  much  as  about  106  ions  are  necessary,  and  even  more  in  certain  specific  cases  like  the  investigation  of  the 
long-time  effects. 

3.RESULTS  AND  DISCUSSION 

3.1.  Recombination  in  ensembles  of  clusters  of  nonseparable  ion  pairs  - short  MFP  case. 

There  are  numerous  examples  of  multiple  ion-pair  systems.  The  ionic  solutions  which  undergo  some  form  of  ionic 
reaction,  ion  pairing  or  nucleation  make  an  obvious  example.  Typical  concentrations  of  ions,  from  10-4  M to  10  M,  cause 
any  reduction  of  the  system  to  a collection  of  isolated  ion  pairs  to  be  rather  far  from  physical  reality.  Another  example  of  the 
ensembles  of  nonseparable  ion  pairs  is  provided  by  the  ionizing  radiation.  A high-energy  particle  when  passing  through 
condensed  medium  leaves  in  its  wake  clusters  of  ion  pairs  - 2,  3 ,4  or  even  more  ion  pairs.  In  the  case  of  heavy  particles, 
like  the  a particles,  we  observe  the  columns  of  ionization.  The  ion  pairs  can  not  be  separated  and  treated  as  the  isolated 
entities.  In  the  study  of  the  chemical  effects  of  high-energy  radiation  the  consideration  of  the  nonhomogen eous  spatial 
distribution  of  the  chemically  active  transients  in  the  track  plays  an  important  role  and  considerable  attention  has  been  given 
to  th  is  problem,  notably  for  the  case  of  water  and  nonpolar  liquids. 

The  recombination  processes  in  spurs  with  more  than  one  ion  pair  were  studied  by  several  authors.  For  the  hydrocarbon 
systems  characterized  by  short  free  paths  of  charged  species,  Bartczak  and  Hummel 6,7  proposed  a random  flight  simulation. 
They  found  that  the  multipair  effects  may  lead  to  substantial  changes  in  the  recombination  kinetics,  the  ion  escape 
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probability  at  infinite  time  as  well  as  the  singlet/triplet  ratio  among  the  recombination  products  8>  9 . Extensive  studies  of  the 
multipair  effects  in  different  systems  were  performed  by  Green  and  coworkers  (see  10,11  and  references  therein). 


>■< 


H 

I — ( 

►J 

I — < 

PQ 

< 

ffl 

O 

Pi 


0- 


$ 

t 1 

g 

2 


CO 


'Z 

O 


-5  -4  -3  -2  -10  12 

lOg(x) 


Fig.  1.  The  ion  survival  probability 
as  a function  of  the  logarithm  of 
time  for  single-pair  (full  curve)  and 
for  multi-pair  systems.  The 
calculations  correspond  to  the  short 
MFP  or  diffusion  model.  The 
broken  curves  correspond  to 
clusters  initially  composed  of  2,  3, 
4,  5 and  10  ion  pairs,  respectively. 
The  initial  distance  between  cations 
of  different  pairs  is  R++= 0.5  nm, 
the  initial  distance  between  cation 
and  anion  of  the  same  pair  is:  R+_ 
=6.0  nm  (upper  figure),  R+.  =8.5 
nm  (middle  part)  and  R+.  =12.3  nm 
(lower  figure).  The  diffusion 
coefficients  for  cations  and  for 
anions  are  D+  = D.=  1.26  X 10'9 
m2/s,  the  Onsager  radius  is  Rc= 28.5 
nm,  the  reaction  radius  #=1.0  nm. 
The  time  is  expressed  as  a 
dimensionless  parameter  T =Dt/Rc. 


The  simulation  of  the  ion  recombination  processes  has  been  performed  for  systems  of  isolated  clusters  of  ions, 
containing  1,2, 3, 4, 5 and  10  nonseparable  ion  pairs 6’ 7 (later  it  was  extended  to  large  systems,  such  as  ionization  tracks  13). 
An  example  of  the  results  of  the  simulation,  in  the  form  of  the  time  evolution  of  the  ion  survival  probability,  is  given  in 
Figure  1.  The  calculations  allow  us  to  compare  the  recombination  kinetics  calculated  for  the  single  pairs  (the  solution  of  the 
Debye-Smoluchowski  equation)  with  those  calculated  for  the  multiple-pair  clusters.  The  initial  distance  between  the  cation 
and  the  anion  of  a pair  is  assumed  as  R+=  6.0  nm,  8.5  nm  and  12.3  nm  for  the  upper,  middle  and  lower  part  of  the  figure, 
respectively.  In  the  case  of  the  clusters  composed  of  more  than  one  ion  pair,  the  cations  are  placed  on  a straight  line, 
R++= 0.5  nm  apart,  and  the  anion  from  a given  pair  is  selected  in  random  direction  around  its  geminate  cation.  We  performed 
the  test  calculations  for  different  geometrical  arrangements  of  the  cations  and  the  results  were  very  similar.  The  calculations 
correspond  to  the  ions  in  a nonpolar  medium,  with  the  dielectric  permittivity  e=2  and  the  Onsager  radius  Rc= 28.5  nm. 
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The  figure  proves  that  the  multiple-pair  effects  in  the  diffusion  equation  are  indeed  significant.  The  full  line  on  the 
right  hand  side  of  the  series  of  the  kinetic  curves  represents  the  single-pair  solution  of  the  DS  equation.  One  can  see  that 
even  the  recombination  of  ions  in  two-pair  clusters  shows  some  deviations  from  the  single-pair  model.  When  we  refer  to  the 
kinetics  obtained  for  the  10-pair  clusters,  the  kinetics  of  the  process  is  shifted  by  a time  decade  (the  abscissa  is  expressed  in 
the  logarithmic  scale)  with  respect  to  the  DS  kinetics. 

Figure  2 shows  the  dependence  of  the  ion  escape  probability  on  the  initial  cation  - anion  distance  for  single  pairs  (lull 
line  calculated  from  the  Onsager  formula)  and  for  the  clusters  of  2,  3,  4,  5 and  10  ion  pairs.  The  abscissa  is  defined  as  the 
inverse  distance  1 /R+_ , the  scale  of  the  ordinate  is  logarithmic.  In  this  coordinate  system  the  Onsager  dependence  is  linear. 
The  dashed  and  dotted  curves  represent  the  best  fitting  to  the  results  for  2, 3,4,5  pairs.  For  these  clusters,  the  curvature  of  the 
fitted  curves  is  very  slight.  The  fitting  for  the  10-pair  clusters  shows  more  significant  deviation  from  linearity.  The  ion 
escape  probability  for  the  clusters  of  more  than  one  ion  pair  is  lower  than  the  Onsager  probability.  It  is  important  to 
underline  that  this  rule  is  true  only  for  short  distance  between  the  cations,  /?++= 0.5  nm,  for  which  the  calculations  of  Fig.  2 
have  been  performed. 


R+_  , nm 


Fig.  2.  The  Pesc  for  the 
diffusion  model  as  a 
function  of  the  initial  cation- 
anion  distance  for  ions  in 
groups  of  2,  3,  4,  5 and  10 
ion  pairs  with  the  separation 
of  #++=0.5  nm  between 
cations  and  for  single  pairs 
(full  line,  calculted  from  the 
Onsager  formula).  The 
parameters  are:  the  reaction 
radius  #=1.0  nm,  #c= 28.5 
nm.  The  distance  #+.  is 
given  as  the  upper  abscissa, 
1/R+.  is  given  as  the  lower 
abscissa. 
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Interesting  example  which  shows  the  potential  of  the  simulation  method  in  analysis  of  ion  diffusion  process  is  given 
below.  Two  important  scaling  properties  of  the  DS  equation  can  be  mentioned:  i)  the  time  evolution  of  the  ion  density  and 
its  spatial  integral,  i.e.  the  ion  survival  probability,  depends  on  the  sum  of  the  diffusion  coefficients  D+  + D.  but  does  not 
depend  on  the  individual  coefficients;  ii)  the  probability  Pesc  that  an  ion  escapes  recombination  after  infinite  time  does  not 
depend  on  the  diffusion  coefficients  of  the  cations  and  anions  at  all. 

Using  the  computer  simulation  method  we  calculated  the  ion  survival  and  the  ion  escape  probability  as  a function  of 
time  in  the  systems  composed  of  independent  clusters  of  pairs,  from  IV=1  to  10  pairs.  The  sum  of  the  diffusion  coefficients 
is,  for  all  the  simulations,  D+  + D = 2D  = 2.528  X 10'9  m2/s.  The  calculation  of  the  ion  escape  probability  as  a function  of  N 
was  performed  for  three  cases:  for  immobile  anions  D=  0 and  mobile  cations,  for  the  same  mobilities  of  cations  and  anions 
D.=D_=D  and  for  immobile  cations  D+=  0.  The  recombination  kinetics  is  different  for  these  cases,  in  contrast  with  the  first 
scaling  property  of  the  DS  equation.  We  observe  quite  significant  difference  between  the  extremal  cases  of  the  escape 
probability  Pesc:  D+= 0,  and  D=  0.  The  ions  in  the  D+=2D  case  recombine  faster  at  the  initial  phase  of  the  process  but  the  ion 
decay  slows  down  at  long  times  and  its  asymptotics,  i.e.  the  escape  probability,  is  higher  than  the  D+=  0 result. 
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3.3.  Kinetics  and  rate  constant  of  bulk  recombination  of  ions  in  solutions 

It  is  assumed  (ref.1,  p.59)  that  for  high  ionic  concentrations  the  rate  of  recombination  decreases  below  the  Debye- 
Smoluchowski  limiting  rate 


kDS  = 4nRcD/{l  - exp(-RJR)  ~ 4tcR,D  (7) 

towards  the  rate  typical  for  neutral  molecules: 

ks  = 4nRD.  (8) 

The  conjecture  is  reasoned  out  as  follows.  In  the  system  of  many  ionic  charges,  the  charges  of  an  individual  pair  are  in 
a sense  neutralized.  When  two  ions  are  alone  in  space  the  Coulomb  field  dominates  the  diffusion  processes  but  in  the  system 
of  a large  number  of  charges  the  ions  should  behave  like  almost  neutral  molecules  and  their  effective  radius  of  influence  is 
reduced  from  a large  Onsager  radius  Rc  to  a much  smaller  radius  comparable  with  the  encounter  radius  R.  In  the  case  of 
hydrocarbons  we  can  expect  quite  substantial  reduction  of  the  rate  constant  at  high  concentrations  (by  a factor  of  100).  For 
polar  solutions  the  reduction  of  the  rate  constant  is  not  so  spectacular,  but  even  in  the  case  of  aqueous  solutions  we  can 
expect  the  reduction  by  a factor  of  2-3 . 

The  recombination  kinetics  described  so  far  is  calculated  for  the  multiple-pair  but  finite  systems.  In  order  to  simulate 
infinite  systems  the  periodic  boundary  conditions  have  been  introduced.  We  place  a certain  number  of  ions  in  the  simulation 
box.  If  we  assume  a given  size  of  the  box  then  the  number  of  ions  is  calculated  from  the  concentration.  The  central  box  is 
then  repeated  periodically  in  all  three  directions:  x,  y,  z.  The  typical  size  of  the  box  differs  from  10  nm  to  20  nm  and  usually 
contains  few  hundred  independent  ions.  The  infinite  size  of  the  system  causes  also  another  problem.  The  Coulomb  forces 
are  of  infinite  range  and  the  ions  from  even  very  far  boxes  interact  with  the  ions  from  the  central  box.  Therefore  the  field 
acting  on  a given  ion  has  to  be  calculated  by  summing  up  the  contributions  from  infinite  number  of  ions.  The  Ewald  method 
has  been  used  to  calculate  these  forces  in  our  simulations.  The  method  is  explained  in  the  Allen  and  Tildesley  book  12. 


Fig.  3.  The  ion  survival 
probability  as  a function  of 
the  logarithm  of  the 
dimensionless  time  x 
obtained  by  the  computer 
simulation  of  the  diffusion 
and  recombination  of  ions  in 
solutions  with  the  dielectric 
constant  e=80  and  for  the 
ion  concentrations  104  M, 
10'3  M,  10'2  M,  10"1  M,  1 M. 


Figure  3 shows  the  evolution  of  the  relative  concentration  of  ions  for  the  ion  recombination  in  ionic  solutions,  for  a 
series  of  initial  concentrations:  104  M,  10'3  M,  10'2  M,  10"1  M,  1 M.  The  simulation  was  performed  for  the  dielectric 
constant  e=80,  which  corresponds  to  an  aqueous  solution.  The  figure  shows  the  recombination  kinetics  as  a function  of  the 
dimensionless  time  x,  which  is  defined  as  z=Dt/Rc2  and,  for  the  assumed  parameters,  we  have  the  relation  between  the  real 
time  and  the  dimensionless  time  that  is  given  in  the  figure.  The  time  axis  is  expressed  as  the  logarithmic  scale.  The  time 
scale  spans  over  eight  time  decades. 
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Fig.  4.  The  results  of  the  rate 
coefficient  calculations  for 
the  ion  concentration  C=1  M 
in  an  aqueous  solution 
(e=80). 

The  circles  and  crosses 
show  the  rate  coefficient  of 
recombination  as  a function 
of  the  logarithm  of 
(dimensionless)  time 

calculated  by  direct 
differentiation  of  the 
simulated  kinetics  of  Fig.  3. 
The  circles  correspond  to  the 
Boltzmann  initial  distribution 
of  cation  - anion  distances, 
the  crosses  correspond  to  a 
random  configuration  of  ions. 
The  full  and  dashed  curve 
show  the  time  dependence 
rate  that  results  from  the 
Debye-Smoluchowski  model 
for  the  Boltmann  and  random 
condition,  respectively,  and 
the  dotted  horizontal  line  is 
the  stationary  rate  constant. 


The  reaction  rate  constant  has  been  calculated  by  direct  differentiation  of  the  ion  survival  probability  P(r)  over  the  time 
x . The  results  of  the  rate  constant  calculation  are  given  in  Figure  4.  These  results  show  the  time  dependence  of  the  rate 
constant  The  cricles  show  the  simulated  rate  coefficient  for  the  Boltzmann  initial  distribution  of  the  cation  - anion  distances 
and  the  crosses  for  the  random  configuration  of  ions.  The  Debye-Smoluchowski  rate  constant  (calculated  from  the  full  time 
- dependent  formula)  is  shown  also  in  the  figure  as  the  full  line  for  the  Boltzmann  condition  and  the  dashed  curve  for  the 
random  condition.  The  asymptotic  value  of  the  rate  is  kDS  = 1.614  x 1010 1/mole/s,  more  than  two  times  higher  than  the  rate 
for  neutral  molecules  ks  = 4kRD  = 0.747  x 1010  1/mole/s  (the  encounter  reaction  radius  has  been  assumed  as  0.4  nm).  The 
simulated  rate  reaches  the  steady  state  value  in  the  last  decade.  The  rate  is  2.05  X 1010 1/mole/s. 

Figures  5 and  6 summarize  the  results  of  the  rate  constant  calculations.  In  Figure  5 the  rate  constants  are  shown  as  a 
function  of  the  initial  concentration  of  ions  in  the  solution,  C(0).  Two  sets  of  the  rate  constants  are  given:  for  e=10  and 
e=80.  The  corresponding  horizontal  lines  show  the  values  of  the  steady-state  Debye-Smoluchowski  rate  constant.  In  Figure 
6 the  rate  constants  are  given  as  function  of  the  dielectric  permittivity  of  the  solution  e.  The  simulated  results  (squares)  for 
C(0)=0.01  M,  for  which  the  interionic  distance  is  on  average  about  4.4  nm,  are  compared  with  the  theoretical  Debye- 
Smoluchowski  rates  (full  curve)  and  with  the  rate  for  the  neutral  molecules  (horizontal  dotted  line). 

Commenting  upon  the  results  of  these  calculations,  we  found  rather  interesting  and  unexpected  behaviour  of  the 
diffusion-controlled  recombination  processes  in  concentrated  ionic  systems.  According  to  qualitative  conjectures,  it  was 
expected  that  the  rate  constant  will  decrease  with  the  increase  of  the  ion  concentration,  from  the  Debye-Smoluchowski 
limiting  rate  for  the  isolated  pairs  to  the  rate  constant  ks=4nRD  typical  for  neutral  molecules  (represented  in  Figures  5 and  6 
by  the  horizontal  dashed  line).  However,  the  picture  obtained  from  the  simulations  is  different.  Roughly,  we  can  summarize 
it  as  follows:  the  ions  of  the  reacting  ion  pair  are  not  influenced  by  all  the  other  charges  of  the  system  and  their  reaction  rate 
does  not  depend  on  the  overall  density  of  charge.  The  discussion  of  these  results  from  the  point  of  view  of  experimental 
measurements  of  the  rate  constant  of  ionic  reactions  is  given  in  17, 19.  Marciniak  and  Rozwadowski  18  investigated  the 
reaction  rate  constants  for  the  recombination  of  the  positively  charged  organic  ions  (derivatives  of  pyridine)  with  a series  of 
singly  and  multiply  charged  ions  with  the  concentration  of  between  10'3  and  10'2  mole/1.  They  obtained  the  rate  constants 
about  1.3  x 1010 1/mol/s  whereas  the  rate  constants  for  similar  reactions  but  with  neutral  molecules  were  about  4 X 109 
1/mol/s,  more  than  3 times  lower.  These  results  are  in  very  good  agreement  with  our  simulation  results. 
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Fig  5.  The  compilation  of  the  rate 
constants  from  the  simulations  of 
the  homogeneous  recombination. 
The  steady  state  rate  constant  is 
given  as  a function  of  the  initial  ion 
concentration  in  the  solution.  The 
squares  show  the  results  for  slightly 
polar  solutions  (s=10),  the 
diamonds  correspond  to  aqueous 
solutions  (e=80).  The  horizontal 
dotted  lines  in  the  vicinity  of  the 
simulated  rate  constants  show  the 
values  of  the  theoretical,  single-pair 
rate  kDS.  The  dashed  line  at  the 
bottom  of  the  figure  shows  the 
value  of  the  rate  constant  for  the 
neutral  particles  reacting  at 
encounter. 


Fig  6.  The  rate  constant  of  the  ion 
recombination  in  the  solutions  with 
the  initial  concentration  C(0)=10'2 
M calculated  as  a function  of  the 
dielectric  permittivity  of  the 
solution.  The  results  of  the 
simulation  are  shown  as  the 
squares,  the  single-pair  Debye- 
Smoluchowski  rate  constant  is 
shown  as  the  full  curve,  the  rate 
constant  for  neutral  molecules  with 
R=0.4  nm  is  represented  by  the 
horizontal  dotted  line. 


3.4.  The  charge  recombination  in  the  systems  characterized  by  anisotropy  of  diffusion  coefficient 

Studies  of  the  migration  of  the  charged  species  induced  by  high-energy  radiation  in  various  solid  systems  have  been 
carried  out  recently  15.  In  some  cases  the  kinetics  of  the  charge  recombination  is  complicated  by  the  feet  that  the  mobility  of 
the  charged  species  differs  in  different  directions.  For  example,  in  crystalline  polyethylene  a large  anisotropy  has  been 
observed,  the  mobility  having  been  found  to  be  at  least  a factor  of  100  larger  in  the  direction  of  the  chain  than  perpendicular 
to  it  l6.  In  various  "conducting"  polymers  the  charge  migration  takes  place  along  the  chain  and  in  phtalocyanine  and 
porphyrin  aggregates  along  the  columnar  stacks. 
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For  the  computer  experiments  discussed  in  this  chapter  14  we  assume  that  the  medium  is  characterized  by  the  diffusion 
coefficient  D,  in  the  z-direction  and  the  diffusion  coefficient  Dx=Dy=D  in  the  direction  perpendicular  to  die  z-axis.  The 
degree  of  anisotropy  is  represented  by  the  coefficient  S=D/D.  The  diffusion  tensor  is  thus  assumed  in  the  diagonal  form, 
without  any  off-diagonal  components.  The  diagonal  elements  Dz  and  Dx=Dy  can  differ  as  much  as  by  factor  of  10  in  our 
calculations,  although  in  molecular  crystals  the  ratio  <5  does  not  exceed  the  value  of  2-3.  In  the  present  calculations  we 
assume  <5  as  a free  parameter  and  we  perform  the  simulations  for  selected  values  of  <5  from  the  range  1-10.  When  the 
calculations  have  to  be  compared  with  the  isotropic  case,  we  calculate  the  effective  diffusion  coefficient,  defined  as  the 
cube  root  of  the  determinant  of  the  diffusion  tensor,  Deg=D  (S)1'3 . 

Figure  7 shows  the  kinetics  of  the  ion  recombination.  The  probability  of  the  ion  survival  against  recombination  up  to 
time  t is  plotted  as  a function  of  the  logarithm  of  the  dimensionless  time  z=Dt/R2c , D=Deff+  + Def.,  and  Deff  is  the  cube  root 
of  the  determinant  of  the  respective  diffusion  tensor.  The  kinetics  in  Fig.  8 were  calculated  for  relatively  high  anisotropy 
coefficient  <5=10.  Different  curves  in  the  figure  correspond  to  different  angles  Q between  the  z-axis  ( i.e . the  direction  of  the 
fast  movement)  and  the  axis  defined  by  the  cation  and  anion  of  the  pair.  The  calculations  were  performed  for  the  angles:  0° 
(the  pair  axis  aligned  with  the  z-axis),  30°,  60°,  90°  (the  pair  axis  perpendicular  to  the  direction  of  the  fast  movement).  These 
calculations  are  represented  by  the  broken  or  dotted  curves.  The  thick  full  line  gives  the  results  of  the  calculations  for  the 
averaged  kinetics:  (he  negative  ion  was  produced  on  a sphere  in  a random  direction  around  the  positive  ion.  In  all  the  cases 
the  initial  cation  - anion  distance  was  fixed:  the  initial  distance  was  R+=  7.7  nm.  Finally,  the  thin  continuous  curve  shows 
the  corresponding  kinetics  in  isotropic  media.  The  diffusion  coefficients  are:  for  cations  D+Z-SD+,  D+X=Z)+,  Dry=  D+-2.521 
X 10"10  m2/s  , for  anions  D.:=SD_,  D.X=D_,  D.y=D  =2.521  X 10'8  m2/s.  The  medium  parameters  correspond  to  a nonpolar 
medium  with  the  On sager  radius  Rc=  30  nm. 

The  recombination  kinetics  for  different  angles  6 are  significantly  different  one  from  another  even  for  low  anisotropy. 
For  the  anisotropy  coefficient  <5=10  the  curves  for  6 =0°  and  #=90°  are  shifted  by  more  than  a time  decade.  The 
recombination  at  low  angles  6 is  the  fastest  and  the  escape  probability  is  the  lowest.  In  the  contrary,  for  high  0 the 
recombination  process  is  slower  and  the  escape  probability  very  high.  For  <5=2  the  escape  probability  at  8=90°  is  almost  4 
times  higher  than  the  probability  at  £H)°,  for  <5=10  the  ratio  of  the  escape  probabilities  can  be  as  high  as  30  - 40.  The  angle 
averaged  kinetics  of  recombination  has  no  longer  typical  shape  of  the  diffusion-controlled  kinetics  but  extends  over  large 
number  of  time  decades  and  becomes  similar  to  the  dispersion-type  kinetics  with  time-dependent  rate  constant. 


Fig.  7.  The  ion  survival 
probability  as  a function  of 
logarithm  of  (dimensionless) 
time  for  the  systems  with  axial 
anisotropy.  The  anisotropy 
coefficient  <5=10.  The  series  of 
broken  curves  correspond  to 
the  following  angles  6 between 
the  z-axis  of  the  coordinate 
system  and  the  position  vector 
of  the  anion  relative  to  the 
position  of  the  cation  of  a pair  : 
0=0°,  30°,  60°  and  90°.  The 
thick  continuous  line  shows  the 
kinetics  averaged  over  the 
angle  6,  the  thin  line  is  the 
isotropic  kinetics. 
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3.5.  Kinetics  and  escape  probability  for  different  lengths  of  MFP  - recombination  of  isolated  ion  pairs. 

Electron  mobility  //  in  chain  or  cyclic  hydrocarbons  is  low:  in  n-pentane  and  n-hexane  at  room  temperatures  is  1.4  x 
10'5  m2/Vs  and  0.9  x 10"5  m2/Vs,  respectively,  in  cyclohexane  3.5  x 10'5  m2/Vs  (see  for  instance 5 ).  There  is  however  a class 
of  liquids  where  an  excess  electron  has  much  higher  mobility:  tetramethylsilane  (jli=1  x 10'2  m2/Vs  at  296  K),  methane 
((p=3.73  x 10'2  m2/Vs  at  120  K),  argon  ((|X=4.9  x 10'2  m2/Vs  at  87  K),  etc.  Using  Langevin  relation  (see  ref. 5)  between  the 
electron  MFP  X and  the  mobility  // : ji=X  (e  /3kBT)  (where  is  the  root  mean  square  velocity  of  thermalized  electron 

calculated  from  the  Maxwellian  distribution)  we  obtain  the  MFP  for  electron  in  n-pentane  or  n-hexane  of  the  order  of  1 x 
10‘2  nm  but  in  tetramethylsilane  the  MFP  is  6.05  nm,  in  methane  22.5  nm  and  in  argon  29.6  nm.  The  latter  MFP's  are  similar 
to  the  Onsager  distance  in  these  media. 

In  the  case  of  a large  MFP  the  electron  may  pass,  keeping  to  an  orbital  trajectory,  at  a very  short  distance  from  a 
positive  ion,  without  recombining,  and  move  away  from  the  cation  again.  Experimental  evidence  for  this  phenomenon  is 
given  by  the  decrease  of  the  recombination  rate  kr  for  high  mobility  systems,  e.g.  liquid  methane.  Also  the  increase  of  the 
yield  of  escaped  charges  in  such  systems  as  liquid  methane,  neopentane  and  tetramethylsilane  in  comparison  with  low 
mobility  liquids  may  be  an  indication  of  this  phenomenon.  An  excellent  review  of  the  experimental  and  theoretical  results 
on  ion  recombination  in  the  systems  characterized  by  long  mean  free  path  of  electrons  is  given  by  Schmidt  5 . The 
calculations  (see  refs. s' 20,  for  a source  of  references)  have  proven  that  when  the  MFP  of  electrons  is  increasing  towards  the 
values  of  die  order  of  the  Onsager  length  Rc  in  a given  medium,  the  theory  of  transport  and  recombination  of  charges  that  is 
based  on  the  DS  equation  does  not  describe  properly  the  recombination  kinetics,  ion  escape  probabilities  nor  rate  constant  of 
recombination.  In  the  case  of  the  geminate  recombination,  the  escape  probability  at  infinite  time  becomes  higher  than  the 
probability  that  results  from  the  Onsager  formula,  much  higher  for  MFP  of  the  order  of  Rc.  The  rate  constant  for  electron  - 
cation  recombination  in  high-mobility  liquids  kr  is  lower  than  the  DS  rate  kDS=4nDRc.  There  is  no  explicit  calculation  of  the 
electron  - cation  geminate  recombination  process  in  the  system  of  more  than  one  e‘-cation  pair.  This  is  particularly 
important  in  the  experiments  where  the  electron  - cation  pairs  were  produced  by  ionizing  radiation  and  a large  part  of  the 
ionization  appears  in  the  form  of  clusters  of  nonseparable  ion  pairs. 

In  order  to  compare  the  results  for  long  MFP  with  the  calculations  performed  for  short  MFP  and  low  mobility 
electrons,  the  first  series  of  the  simulations  has  been  performed  for  single  pairs  electron  - cation.  The  results  of  the 
calculations  of  the  kinetics  of  recombination  are  presented  in  Figure  8.  The  kinetic  curves  are  plotted  versus  logarithm  of  the 
dimensionless  time  x defined  as  r=Dt/Rc2  where  D-e/j/kBT  is  the  diffusion  coefficient  of  the  electrons  (the  diffusion 
coefficient  of  the  cations  is  negligible  as  compared  to  D). 


Fig.  8.  The  probability  of  electron 
survival  up  to  the  time  t plotted  as 
a ftmction  of  the  logarithm  of 
T=Dt/Rc2.  The  calculations  were 
performed  for  single  ion  pairs.  The 
different  curves  correspond  to 
different  length  of  the  mean  free 
paths  (MFP)  of  electron:  from  X=0 
(the  infinitely  short  MFP  or  the  DS 
case),  for  X=0.05RC,  0. 107?c, 

0.20J?C,  0.35 Ra  0.50 Rc,  0.75 Rc  to 
A,=1.0E„  where  Rc- 28.5  nm  is  the 
Onsager  radius.  The  cation  - 
electron  initial  distance  is  R+. 
=0.45  Rc. 
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The  full  curve  is  calculated  on  a basis  of  the  diffusion  model  and  represents  the  DS  kinetics  of  the  ion  recombination 
for  isolated  ion  pairs  in  a nonpolar  medium  with  Rc= 28.5  nm.  Physically,  this  kinetics  corresponds  to  die  electrons  with 
mobility  of  the  order  of  5 x 10'8  m2/Vs,  i.e.  the  mobility  of  the  classical  ions.  The  remaining  kinetic  curves  correspond  to  the 
mean  free  path  of  the  electrons  of  0.05iJc , 0.10/?„  0.20i?„  0.35 R„  0.5QRC,  0.75 Rc  and  1.00#c  or  to  the  electron  mobilities 
from  2.4  x 10'3  m2/Vs  to  4.7  x 10‘2  mWs  for  0.05 Rc  and  1 ,007?c , respectively.  It  is  worth  to  note  a dramatic  difference 
between  the  DS  kinetics  and  the  kinetics  for  the  MFP  1=1 .007?c=28.5  nm.  The  results  for  die  MFP  1=0.20RC  are  also  very 
much  different  from  the  diffusion-controlled  limiting  case.  The  DS  kinetics  may  be,  with  acceptable  errors,  applied  only  to 
the  case  of  A=0.05/?C=T.43  nm  and  this  seems  to  be  the  limit  of  applicability  of  the  DS  equation. 

Figure  9 shows  the  results  of  the  calculations  of  the  probability  that  an  electron  - cation  pair  ultimately  escapes 
recombination.  The  probability  is  plotted  as  a function  of  the  initial  cation  - electron  separaton  R+.  .The  results  were 
obtained  from  the  long-time  parts  of  the  kinetic  curves  by  fitting  the  kinetics  to  the  function  and  taking  the  limit 
T— *».  A series  of  calculations  were  performed  for  the  following  values  of  the  MFP  1:  0.05 Rc,  0.10/?o  0.20J?C,  035RC, 

0. 50RC,  0.75 R„  1.00/Jo  labeled  in  the  figure.  The  curve  marked  with  "+"  and  the  label  ONS  was  calculated  from  the 
Onsager  formula  for  the  diffusion-controlled  case:  Pesc=  exp(-Rc/R+ _).  Similarly  as  in  the  case  of  the  kinetic  curves  of  Figure 

1,  we  can  see  dramatic  differences  between  the  Onsager  escape  probability  and  the  probability  for  the  MFP  X>0.05RC . The 
increase  of  the  escape  probability  above  the  Onsager  limit  exceeds  an  order  of  magnitude. 


R+.  , nm 


Fig.  9.  The  escape 
probability  as  a function 
of  the  inverse  of  the  initial 
cation  - electron  distance 
R+_.  The  calculations  were 
performed  for  single  ion 
pairs.  The  different  curves 
correspond  to  different 
lengths  of  the  mean  free 
paths  (MFP)  of  electron: 
from  1=0  (the  Onsager 
case,  marked  ONS),  for 
T=0.05.RC,  0.1i?c,  0.2 R„ 
0.35 Rc,  0.5 Ra  0. 75 .ft c, to 
1 .Rc. 


4.  CONCLUSIONS 

The  limits  of  validity  of  the  Debye  - Smoluchowski  equation  applied  to  diffusion-controlled  recombination  of  charges 
were  investigated  by  means  of  the  computer  simulation  method. 

The  scale  of  the  deviations  of  the  kinetics  of  the  recombination  process  in  the  multi-pair  clusters  from  the  kinetics  for 
the  isolated  pairs  was  estimated.  The  deviations  of  the  multi-pair  kinetics  and  escape  probability  from  the  corresponding 
single-pair  results  are  significant.  The  kinetic  curve  for  recombination  in  multi-pair  clusters  may  be  shifted  by  more  than 
time  decade  from  the  DS  kinetics  and  the  escape  probability  can  be  lower  by  an  order  of  magnitude  than  the  Onsager 
probability  for  the  same  set  of  parameters.  However,  in  the  case  of  the  bulk  recombination  of  ions  in  solution,  in  spite  of  the 
significant  difference  of  the  simulated  kinetics  and  the  DS  kinetics,  the  steady-state  rate  constants  of  recombination  are  quite 
close  to  the  k^  rate,  even  for  very  high  concentration  of  ions. 
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Some  differences  between  the  simulated  kinetics  and  the  DS  kinetics  were  also  noted  for  the  systems  with  marked 
anisotropy  of  the  diffusion  tensor. 

For  the  processes  with  long  mean  free  path  of  charges  the  scale  of  deviations  is  much  higher.  If  the  length  of  the 
electron  MFP  exceeds,  say,  0.057?„  the  deviations  of  the  recombination  kinetics  and  the  escape  probability  from  the 
respective  values  obtained  on  a basis  of  the  diffusion  mechanism  are  so  large,  that  neither  the  DS  equation  nor  the  Qnsager 
formula  can  be  applied  to  describe  the  recombination  between  positive  and  negative  charges. 
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